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Abstract —We introduce a new numerical method for the 
computation of the inverse nonlinear Fourier transform and 
compare its computational complexity and accuracy to those of 
other methods available in the literature. For a given accuracy, 
the proposed method requires the lowest number of operations. 
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transform (EFT) algorithm. All the methods are tested for 
the focusing NLSE in the soliton-free case (i.e., when only 
the continuous spectrum of the Zakharov-Shabat operator 
is present), which is of particular interest for information 
transmission based on nonlinear inverse synthesis a. Finally, 
we investigate and compare the accuracy and computational 
complexity of all the methods. 


I. Introduction 


II. Numerical methods 


The propagation of light in optical fibers is governed by 
the nonlinear Schrodinger equation (NLSE), which accounts 
for the interplay of dispersion and nonlinearity. In the lossless 
case, the NLSE admits an analytical solution based on the 
inverse scattering transform (also known as nonlinear Eourier 
transform—NET) |[T|. On this basis, during the ’80s and the 
’90s, research on soliton transmission in optical fibers was 
carried out extensively. However, despite a great effort, soliton 
systems never caught on commercially and were eventually 
abandoned in favor of simpler and more effective transmission 
techniques originally developed for linear channels. Recently, 
an alternative approach to communication over nonlinear 
channels—originally proposed in ID and further generalized 
in lO and ill —is gaining new attention from the optical fiber 
community. Rather than using solitons as information carriers, 
the idea is that of encoding information onto the discrete 
and/or continuous spectrum of the Zakharov-Shabat operator 
(the scattering data), which evolves trivially and linearly along 
the fiber, such that the encoded information can be easily and 
accurately recovered from the optical signal after propagation 
at any distance, without any impact from dispersion and 
nonlinearity. This approach requires the computation of the 
inverse NET at the transmitter to obtain the transmitted time- 
domain signal from the modulated scattering data, and the 
computation of the direct NET at the receiver to extract the 
scattering data from the received time-domain signal. Both 
operations are, in general, quite involved and research on 
efficient numerical methods for their implementation is in 
progress a-in. 

In this work, we focus on the numerical computation of 
the inverse NET. We consider two methods available in the 
literature 11, 0 and propose an alternative method based 
on iterated convolutions evaluated through the fast Fourier 


We consider the focusing NLSE in the normalized form 

1 , ,9 

jU:, + -Utt + U\u\ = Q (1) 

We assume to know the scattering data at a given distance 
z and we focus on the inverse scattering problem (or inverse 
NET) to obtain the corresponding NLSE solution u{t, z)Q This 
problem leads to the Gel’fand-Levitan-Marchenko equation 
(GLME) Oa 


t t 

K{t,y) + F{t + y)-\- J J K{t,r)F*{r + s)F{s + y)dsdr = 0 
— 00—00 

( 2 ) 

where the integral kernel F{y) is a function of the scat¬ 
tering data. Specifically, considering the soliton-free case, 
F{y) = ^ r(X)e~^^'^dX where r(A) is the reflection 

coefficient. The solution of the NLSE is finally obtained as 
u{t) = 21imj^^t- K{t,y). 

Let F'{y) = F{y) for y > 0 and F'{y) = 0 otherwise. The 
GLME can be rewritten in the form of the Marchenko integral 
equations 
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Bi{t,a) = — J F'*{2t—a—(3)B2{t, P)df3 
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B 2 {t,a) = J F'{2t—a—(3)Bi{t, P)df3 — F'{2t—a) 

0 

and the NLSE solution obtained as u(t) = 2i?2(f,0+). In 
the following we assume to know F on a uniform grid in 
the interval [0, 2T] with step 6a = ^ and we analyze three 
different numerical methods to solve the system 0, in order 
to find the solution of ([III for f G [0,r]. Specifically, we 


^Hereafter, the dependence on 2: is inessential and will be omitted. 


consider the uniform grid on [0, T] given hy tm = {m — l)St 
for m = l,..,Nu + 1 with discretization step 6t = 

The complexity of the methods is measured by considering 
the required number of complex products, assuming that the 
FFT of a vector of N elements requires ^ log 2 (-^) complex 
products. 

An efficient numerical way to compute the inverse NFT, 
assuming Np = Nu, and therefore St = is presented in 
0. The authors, using the Nystrom method with a rectangular 
quadrature scheme to approximate the integrals in Q, find the 
solution u{tm) by solving a linear system characterized by a 
Toeplitz matrix of size Np+m. The method of Trench for the 
inversion of non-Hermitian Toeplitz matrices ca is used to 
recursively compute the solution for every value tm- For this 
reason, let us name this method Nystrom-Trench (NT). The 
resulting algorithm computes the solution u{t) in the desired 
interval with complexity^ 

C ~ 117V| (4) 

However, due to the recursive nature of the algorithm, the com¬ 
plexity remains almost unchanged if the solution is computed 
in a lower number of points Nu- 

Another technique to solve system 0 with St = for 
c S N, (i.e., in = Np/c points in the interval [0,T]) is 
presented in Q, where the Nystrom method with composite 
Simpson’s quadrature rule is considered. For a fixed time tm, 
the resulting linear system is solved by means of the conjugate 
gradient method. Therefore, we refer to this scheme as the 
Nystrom conjugate gradient (NCG) method. Since we deal 
only with diagonal and upper-left triangular Hankel matrices, 
the complexity of the method can be reduced by using the 
FFT algorithm to compute the matrix-vector products required 
by the conjugate gradient method. As a result, the number 
of operations required to compute the solution in the desired 
interval iJ^ 

C ~(1 -I- fcmax) ^ 6[c(to - 1) -f 1] log2{2[c(m - 1) -I- 1]}+ 

m—1 

+ NuNp{A + Qk 

max) 

(5) 

where fcmax is the number of iterations of the conjugate 
gradient method, which is assumed to be the same for every 
t for the sake of simplicity. In this case, the solution is 
independently evaluated at each time tm in the given interval, 
such that C depends significantly on the number of points N^. 

Now, let us define a new iterative method to solve the system 
(O for t S [0,T] with discretization step St = c^. The main 
idea is that, by defining the auxiliary functions Bt{t,a) = 
Bi{t,a) for a > 0 and B[{t,a) — 0 otherwise, with i — 
1, 2, the integrals in Q can be seen as convolutions between 
these functions and the kernel F and efficiently computed by 
means of the FFT algorithm. Therefore, considering a proper 

^The computational complexity of the methods presented in (4j and 
is not indicated hy the authors. The values reported here refer to the most 
efficient implementations that we were able to devise. 


Starting value for Bi and iteratively updating the values of 
B 2 and Bi by alternately computing the integrals in (O, we 
define an iterative method that, under certain conditions on F 
and t, converges to the solution of the system. As a starting 
point for Bi at tm, we pick the solution found at the previous 
point tm-i, as we are dealing with continuous functions. The 
resulting algorithm is 


For fc = 0 
B[°\a,tm) = I 
For fc = 1, .., fctnax 
I B^^\a,tm) = 
\ B^^\a,tm) = 


0 if TO = 1 

B^^““\a,tm-i) if TO > 1 

-F{2tm-a) + {F^B'f-^'>){2tm-a) 
-{F**Bf^){2tm-a) 

(6) 


(k) 

where * denotes convolution, H is the fcth estimate of the 
function Bi, and B \ the corresponding auxiliary function. 
We refer to this method as the iterative convolution (IC) 
method. The number of operations required to find the solution 
in the points inside the interval [0,T] is 

Afu-fl 

C — E 3 [ 2 c(to— log 2 [2c{m—l)+l]+2k^^,^NuNp 

m—1 

(7) 

Also in this case, C depends significantly on the number of 
points Nu and can be reduced by evaluating the solution with 
a larger step size. 


III. Numerical results 


With the purpose of comparing the methods presented 
above, we consider the single-pulse signal considered in 0 


u{t) = 


‘iaiya{a — 1) 

(ct - l)2e-2<^at + jy2g2<Tat 


( 8 ) 


where a = + 1- In this case, the scattering data are 

known and the GLME kernel is 0 


F(y) = ave-°^y (9) 

for a > 0 and — 1 < < 1. We consider the interval [0,3] 

and show the results only for a = 12 = 1, as similar results 
are obtained for different values of the parameters. 

In Fig. [T] we compare the analytical solution ® with the 
numerical results obtained with <5^ = 0.01 and St = 0.005. 
At this scale, all the numerical results are practically super¬ 
imposed (and are, therefore, represented by a single curve 
to avoid cluttering the figure) and equal the exact solution. 
Since the new IC method is an iterative method, whose 
starting point at tm depends on the previous solution found at 
tm-i, it is interesting to investigate its convergence properties 
when the solution is evaluated with a different resolution St 
(keeping unchanged the resolution Sa for the evaluation of 
the integrals). Fig. |2| shows the modulus of the eiTor, i.e., the 
difference between and the numerical approximation, as 
a function of t for Sa = 0.01, different iteration numbers fc, 
and different resolutions St- In all the cases, three iterations 
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Fig. 1. Comparison between the exact solution and the solution obtained Fig. 3. RMSE versus number of iterations for the NCG and the IC method 



Fig. 2. IC method: modulus of the en'or for different iterations k, different 
resolutions <5t, and a fixed resolution (5 q: = 0.01. 


are enough to reach the convergence. However, in the case 
5t = 5a, the error after the first iteration is significantly lower 
than the error obtained after the convergence is reached. This 
quite surprising result has been observed also for different 
resolutions 5a and for different pulse shapes (not shown here 
due to a lack of space). However, we have not been able to 
find a good explanation for this behavior, which is still under 
investigation. Given its particular properties, in the rest of the 
paper we will refer to the IC method with fcmax = 1 and 
5t = 5a as the ICl method. 

Fig. [2 compares the convergence properties of the IC 
method (for 5t = 5a = 0.01) to those of the NCG method 
(for 5a = 0.01 and any 5t, as they are independent of 5t), 
showing how the root-mean-square error (RMSE) evaluated 
over the whole interval [0,3] evolves with the number of 
iterations. Thanks to the composite Simpson’s quadrature rule 
employed by the NCG method (which is more accurate than 
the rectangular quadrature rule implicitly employed by the IC 
method when using the EFT to compute the convolutions), the 
final RMSE reached by the NCG method is lower than the 
RMSE reached by the IC method. However, the IC method 


converges faster (after 3 iterations) than the NCG method 
(after 6 iterations). The same fast convergence of the IC 
method is observed also for higher values of 5t (as shown in 
Fig. |2l. However, in the special case 5t = 5a, as observed 
before, the IC method achieves an even lower RMSE at 
the first iteration (ICl method). Such a fast convergence is 
obtained thanks to the good starting value for Bi indicated 
in (HJ. A convergence of the iterative method is observed in 
all the considered cases. This is, however, not guaranteed in 
general, as the convergence depends on specific conditions on 
F and t that are currently being investigated. 

For all the methods, both the accuracy and the complexity 
depend on the resolution 5a employed for the quadrature rules. 
Fig. m shows the RMSE as a function of 5a for all the con¬ 
sidered methods. In the NT and ICl methods, the resolution 
5t is set to (5a/2 and 5a by construction, respectively. On the 
other hand, in the NCG and IC methods, 5t can be selected 
as an integer multiple of 5a/2 without affecting the accuracy. 
The RMSE values obtained with the NT and IC methods are 
practically the same, as they are both based on a rectangular 
quadrature rule for the computation of the integrals in ([2l. 
The NCG method, thanks to the more accurate composite 
Simpson’s quadrature rule, achieves a lower RMSE for the 
same resolution 5a- Finally, the ICl method achieves an even 
lower RMSE. 

Fig. a can be used to determine, for each method, what is 
















Fig. 5. Number of operations required by the four methods for a fixed RMSE 
of 2 X 10“^ and a variable desired resolution 5t. 

the resolution Sa (and, correspondingly, the number of points 
Np = 2Tj6a) required for the computation of the integrals 
in order to find the solution u{t) with a given accuracy. This, 
in turn, will affect the complexity of the method according 
to (|4|i, Q, or (|7]i. Another relevant parameter that affects the 
computational complexity of the NCG and IC methods (but not 
that of the NT and ICl methods) is the resolution St = cSal2 
with which the solution u{t) is needed (or, equivalently, the 
number of points Nu = Np/c in which the solution must 
be computed). Fig. |5] compares the computational complexity 
required by the various methods, computed according to (|4]l, 
dSll, and (|7]i, to obtain the solution u{t) with a desired RMSE 
of 2 X 10“^ (about in Fig. |4li, as a function of the desired 
resolution St- According to Fig. |4l the required resolution 
is 6a — 0.0086 for the NT and IC methods, Sa — 0.014 
for the NCG method, and Sa — 0.033 for the ICl method. 
The number of iterations required to reach the convergence 
is fcmax = 6 for the NCG method and fcmax = 3 for the 
IC method. By construction, the minimum resolution St is 
(5a/2 for the NT, NCG, and IC methods, and Sa for the 
ICl method. However, if desired, the resolution St can be 
increased by computing the solution u{t) in a lower number 
of points Nu in the given interval [0, T]. This will significantly 
reduce the computational complexity of the NCG and IC 
methods, but will leave the complexity of the NT and ICl 
methods almost unaffectedly For any resolution St, the NCG 
and IC methods have a similar complexity. In particular, the 
latter is slightly more convenient, though the former could be 
probably improved by using a more accurate quadrature rule or 
employing preconditioning. Both methods are more complex 
than the NT method when employed at full resolution, but 
become less complex when the solution is required in a lower 
number of points (the IC method for St > 0.05 and the 
NCG method for St > 0.075). Finally, the ICl method is 
significantly less complex than the NT method (by about one 
order of magnitude) and less complex than the NCG and IC 
methods even at high resolutions St- However, the behavior 

^Even if the solution is needed in a lower number of points Nu < Np, 
the NT and ICl algorithms must be anyway executed for Np points. In 
the NT case, however, some operations can be saved, slightly reducing the 
computational complexity with respect to what is reported in E). 


of the ICl method is still under investigation and its superior 
performance, though observed also in other cases (not reported 
in this paper for a lack of space), has not been fully explained 
nor demonstrated for a generic signal. 

IV. Conclusions 

A new numerical method for the computation of the inverse 
nonlinear Fourier transform has been proposed. The solution 
is obtained after iterated convolutions with the GLME kernel, 
which are efficiently computed trough the EFT algorithm. 
The accuracy and computational complexity of the proposed 
method have been investigated and compared to those of two 
other methods available in the literature, both exploiting the 
Nystrom method in combination with either the conjugate 
gradient algorithm 0 or the Trench algorithm i). The ob¬ 
tained results show that, for a desired accuracy, the proposed 
method requires the lowest number of operations. The relation 
between the discretization step with which the integrals of the 
GLME are approximated and the time resolution with which 
the solution is computed (which need not be the same) is 
also discussed for all the methods, and the dependence of 
accuracy and complexity on those two parameters is inves¬ 
tigated. When the two parameters are equal, the proposed 
method shows a very good but peculiar behavior which is still 
under investigation. Further investigation is also required to 
compare the algorithms in different scenarios, such as for the 
defocusing NLSE or in the presence of discrete eigenvalues 
in the nonlinear spectrum of the signal. 
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